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Program summary 



Program Title: Helac-Phegas 
Journal Reference: 
Catalogue identifier: 
Licensing provisions: none 



Program obtainable from: http : //helac-phegas . web . cern . ch/helac-phegas/ 



Distributed format: tar gzip file 
Programming language: FORTRAN 
Computer: All 
Operating system: Linux 

Keywords: Dyson-Schwinger equations, recursive algorithms, automatic evaluation of helicity 
amplitudes and total cross sections. 
PACS: 12.38.Bx, 13.85.Dz, 13.85.Lg 
Classification: 11.1, 11.2 

External routines/libraries: Optionally Les Houches Accord (LHA) PDF Interface library. 
Nature of problem: One of the most striking features of final states in current and future colliders 
is the large number of events with several jets. Being able to predict their features is essential. 
To achieve this, the calculations need to describe as accurately as possible the full matrix ele- 
ments for the underlying hard processes. Even at leading order, perturbation theory based on 
Feynman graphs runs into computational problems, since the number of graphs contributing to 
the amplitude grows as n\. 

Solution method: Recursive algorithms based on Dyson-Schwinger equations have been developed 
recently in order to overcome the computational obstacles. The calculation of the amplitude, us- 
ing Dyson-Schwinger recursive equations, results in a computational cost growing asymptotically 
as 3", where n is the number of particles involved in the process. Off-shell subamplitudes are 
introduced, for which a recursion relation has been obtained allowing to express an n— particle 
amplitude in terms of subamplitudes, with 1—, 2—, ... up to (n — 1) particles. The color con- 
nection representation is used in order to treat amplitudes involving colored particles. In the 
present version Helac-Phegas can be used to efficiently obtain helicity amplitudes, total cross 
sections, parton-level event samples in LHA format, for arbitrary multiparticle processes in the 
Standard Model in leptonic, pp and pp collisions. 
Does the new version supersede the previous version: yes, partly 
Reasons for the new version: substantial improvements, major functionality upgrade 
Summary of revisions: color connection representation, efficient integration over PDF via the 
Parni algorithm, interface to LHAPDF, parton level events generated in the most recent LHA 
format, k± reweighting for Parton Shower matching, numerical predictions for amplitudes for 
arbitrary processes for phase space points provided by the user, new user interface and the pos- 
sibility to run over computer clusters. 

Running time: Depending on the process studied. Usually from seconds to hours. 
References: 

[1] A. Kanaki, C. G. Papadopoulos, Comput. Phys. Commun. 132 (2000) 306 
[2] C. G. Papadopoulos, Comput. Phys. Commun. 137 (2001) 247 
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1 Introduction 



Current and forthcoming high energy colliders allow the study of final states with many 
hard and well separated jets as well as charged leptons. These multiparton final states 
carry the signature of the decay of massive particles, such as for example W and Z gauge 
bosons, top quarks, Higgs bosons, and possibly supersymmetric particles, as well as the 
effect of new interactions. Therefore the accurate description of such processes is of 
fundamental importance for the study of the properties of the Standard Model particles 
as well as for the discovery of new physics. 

Multiparton processes can be described within standard perturbation theory. At the 
leading order (LO) one needs to calculate matrix elements and efficiently integrate them 
over the available phase space. Although the procedure is well defined and straightforward 
in principle, several challenges are present in practice. First of all, if the matrix elements 
are calculated using the standard Feynman diagrams, the number of which grows facto- 
rially with the number of partons in the process, the set up of an efficient calculational 
framework in the multiparton case becomes impossible. The solution to this problem is 
by now well known and relies on the use of recursive equations instead of Feynman graph 
expansion [1-14]. Helac [9,10] has been the first implementation of the algorithm based 
on Dyson-Schwinger equations. It is able to calculate matrix elements for arbitrary scat- 
tering processes. A fundamental issue is the consistent inclusion of all color configurations 
within QCD for processes with colored partons. The color connection representation is 
used in Helac, and the programme is actually able to calculate all possible color connec- 
tions for arbitrary processes. For amplitudes with many colored particles, a Monte-Carlo 
treatment over color should also be used. This issue has also been addressed so far in 
references [8,11]. Moreover, a new version of the Helac algorithm [12], based on explicit 
color configurations, has been proved very efficient in dealing with processes with many 
colored partons. The next challenge is the phase space integration. For the phase-space 
generation, Phegas [15] has been the first implementation of a complete automated algo- 
rithm of multichannel phase space mappings for arbitrary number of external particles. It 
uses information generated by Helac and automatically performs a multichannel phase 
space generation, utilizing all 'scalarized' Feynman graphs, for any given process, as phase 
space mappings. Finally, since colliding particles like protons are not elementary one has 
to take into account a number of elementary parton level subprocesses. 

It should be mentioned that a number of computer programs, like for example Alpgex 
[17,18], [19], MadEvent [20-22] Sherpa [23,24], Whizard [25] are available, which 
implement some of these techniques and can be used to obtain leading order predictions 
for processes with a rather high number of final state particles. 

In the present paper we will summarize the new developments in the Helac-Phegas 
which fall into two categories: 

1. Functionality: 
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(a) Full implementation of QCD in the color connection representation 

(b) Automatic summation over subprocesses for pp and pp colliders and e'^e~ col- 
liders 

(c) Built-in CTEQ611 PDF support [26,27] with alternative interface to the LHAPDF 
library [28] 

(d) Ut reweighting for matching algorithm to parton shower programs, e.g. Pythia 
[29,30], Herwig [31] or HERWIG++ [32] 

(e) Parton-level event files in the LHA format [33,34] 

(f) Extended set of standard cuts 

2. User interface: 

(a) Program control using a single user input file managed by scripts 

(b) Easy access to numerical amplitude values for a user given phase space point 

(c) Set up for parallel runs on computer clusters 

The present paper is organized as follows. In the next section, after a brief summary of 
the Phegas phase space generator and the Helac algorithm, we describe the treatment 
of the color degrees of freedom, based on the color connection representation and the 
matching to the parton shower algorithms. The rest of the publication describes in detail 
the usage of the program. We summarize with a short outlook. 

2 The Helac-Phegas Algorithm 

The core of Helac-Phegas is made of a phase space generator and a matrix element 
evaluator. The algorithms of both have already been described previously in [15] and [9] 
respectively. Therefore, we will only give here a very short overview of the main ideas. 

2.1 Phegas - Phase Space Generation 

Phegas [15], is a multichannel self optimizing phase space generator. It automatically 
constructs all possible phase-space mappings in order to optimally describe the peaking 
structures of a given scattering process. These phase-space mappings, called also chan- 
nels^ have a one-to-one correspondence to the Feynman diagrams contributing to the 
given process. The information concerning the pole structure of all Feynman diagrams 
is produced by Helac by using the skeleton construction of the solution of the Dyson- 
Schwinger recursive equations. The self-optimization is performed according to the ideas 
presented in [16]. Moreover since in several practical applications the number of phase- 
space channels may become quite large (of the order of many thousands), optimization 
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may be used to significantly reduce this number, based on the contribution of the corre- 
sponding channels to the total variance. At the end a few tens of channels are kept only, 
which is usually sufficient to give an adequate phase space efficiency. Let us stress that 
the selection of these channels is done automatically by the optimization procedure. 

In the case of pp and pp collisions, besides phase space generation for a given center- 
of-mass energy, an integration over the fractions Xi and X2 of the momenta of the initial 
partons, weighted by the parton distribution functions, is necessary. This integration is 
optimized via the Parni algorithm [35]. Its efficiency has been proven to be very good. 

2.2 Helac - Amplitude Computation 

Helac uses an alternative [6-8,11] to the Feynman graph representation of the scattering 
amplitude which is provided by the Dyson-Schwinger approach [9]. 

Dyson-Schwinger equations express recursively the n-point Green's functions in terms 
of the 1—, 2—, . . ., (n — l)-point functions. For instance in QED these equations can be 
written as follows: 




i=l 

+ (^^)^^;V'(W^(A)6(Pl,P2) 

P=Pl+P2 

where 

describes a generic n-point Green's function with respectively one outgoing photon, fermion 
or antifermion leg carrying momentum P. H^,^ stands for the boson propagator and e takes 
into account the sign due to fermion antisymmetrization. 

The computational cost of Helac-Phegas's algorithm grows as ~ 3", which essen- 
tially counts the steps used to solve the recursive equations. Obviously for large n there is 
a tremendous saving of computational time, compared to the n\ growth of the Feynman 
graph approach. 

Helac-Phegas is using a top-down approach in solving the recursive equations. In the 
first phase, a skeleton is constructed, composed by a series of vertices that are participating 
in the given amplitude. All this information is encoded using integer arithmetic. Then this 
information is used in order to actually calculate the amplitude. In this second phase, 
momenta from Phegas are used to compute wave functions. Based on the skeleton 
information the subamplitudes are computed using floating arithmetic. The result is the 
value of the matrix element for each given helicity and color connection configuration. 
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2.3 Color representation 

The treatment of the color degrees of freedom is a very important issue, especially for 
amplitudes involving many colored particles. In the most commonly used approach we 
can express a n— gluon amplitude in terms of the so-called color-ordered amplitudes, using 
the following equation, 

a = ((Ti, o"2, . . . , cr„) is permutation of (1, 2, . . . , ra). The sum is extended over all {n — 1)! 
permutations, where are functions of momenta and helicities of the external particles. 
The important thing is that these Aa color-ordered amplitudes can be calculated using the 
so-called color-ordered Feynman rules. In the framework of the Berends-Giele recursive 
equations, which is the Dyson-Schwinger equations subject to the restriction of ordering 
of the external particles (according to the permutation a), the computational cost grows 
like n^, exhibiting a polynomial growth. When calculating the full color-summed matrix 
element squared, all (n — 1)! terms have to be computed. This approach is some times 
referred as the SU {N^) approach. The inclusion of quark-antiquark pairs is possible, but 
the analogue of Eq.([T]) becomes more complicated. 

The color connection representation is based on a U{Nc) representation of the fields. 
In that case amplitudes involving quarks and gluons can by described in a unified way. 
To this end all indices of the adjoint representation are transformed to indices in the 
fundamental representation by multiplying with t^j, 

Mff;-'- = y . . . t^". . (2) 

JlJ2---Jn Z — d *1J1 injn V / 

a 

When quarks are only involved no action is required. The U{Nc) representation has been 
introduced in Ref. [36]. It has been used in connection with matrix element evaluation in 
Helac-Phegas [9,10]. A Lagrangian based approach has also been used in [37]. 
Now the amplitudes can be written as follows: 

a 

In the above expression we have chosen to fix the order of anti-color indices. Without 
altering the result, the same can be done with the color indices. We have kept the same 
notation for the functions Aa in Eq.([T]) and Eq.([3l). This is indeed the case for the n— gluon 
amplitude but not when quarks are involved. In this representation there is no distinction 
anymore between gluons and quarks. Gluons are represented by a pair of color-anti-color 
lines, where quarks (anti-quarks) are represented by color (anti-color) lines. The number 
of permutations is in principle equal to n^!, where rii is the number of color (anti-color) 
lines. In a process with rig gluons and riq quarks (and anti-quarks), ni = Ug + rig. 

The important thing is again that the A^ function can be calculated using a given set 
of Feynman rules. For instance, the three-gluon vertex can be rewritten as 

'^tAB^CD^EF = ~~:{^adScfSeb — ^af^cb^ed) (4) 
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where on the right hand side only products of 5's appear. Let us introduce a more compact 
notation and associate to each gluon a label {i, ai), which refers to the corresponding color 
index of the previous equation, namely 1 — > A, cti — > S and so on. The use of CTj labels 
will become clear in a moment. With this notation, the first term of the above equation 
is proportional to 

(5) 

and its diagrammatic representation is given by 



The momentum part of the vertex is still the usual one, 

5^12 (Pl - P2f + 923 (P2 - Psf + 931 {P3 " Pi)' 

where gi2 stands for 9n^^2^ ^tc. For the qqg vertex we have 

Xl^AS^CD = -ji^AD^CB - J^^AB^Cd) 



(6) 



(7) 



and if we associate a label (i, 0)((0, ai)) to a quark(antiquark) we have two possible vertices 



(0,a.) 



as well as the one with the 'neutral' gluon, namely. 



(i,0) 



g 



(0,0) 5. 



la- 




(0,CJ.) 



(i,0) 

q 



the last one with an extra factor proportional to i/^/Nl. Finally the four- gluon vertex is 
given by a color factor proportional to 



^10-3 ^3(72 ^2(74 '^4<T1 



(8) 
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and a Lorentz part 



25'125'34 — 913924 — 9u923 



(9) 




(4,a4) 



To make use of the color representation described so far, let us assign a label {i, o'j{i)) to 
each external giuon, (i, 0) to a (incoming) quark and (0, <Ji{i)) to an (incoming) antiquark, 
where i — 1 . . .ni and ai,I — 1 . . . n^! being a permutation of {1 ... n^}. It is clear that 
since all elementary color factors appearing in the color decomposition of the vertices are 
proportional to 5-symbols, the total color factor can only be given by 

= Siaj{l)S2ai{2) ■ ■ ■ Snai{n) (10) 

Moreover the color matrix, defined as 

Cij^ E ^i^l (11) 

colours 

with the summation running over all colors, 1 . . . Nc, has a very simple representation 

Cij = {N^r^"""'^ (12) 

where 1 < m{(7i, uj) < n counts the number of cycles made by the elements of the 
permutations (7/ and Uj. 

The practical implementation of these ideas is rather straightforward. Given the in- 
formation on the external particles contributing to the process, we associate color labels 
of the form (i,cri), depending on their flavor. The Feynman rules to build up the higher 
level subamplitudes, are the ones described above. The result, for each permutation, is 
a new colored ordered ^-function that corresponds to the given color factor J^. The 
computational cost is again given, as in the case of Berends-Giele recursive equations, 
by n^. The total color-summed squared amplitude is obtained after summing over all ni\ 
color connection configurations using the color matrix C. 

In the generation of unweighted events a color connection information is needed for a 
proper matching with parton shower algorithms. For the selection of the color connection 
structure to be assigned to a given event we consider only those v4-functions that corre- 
spond to color connections without a 'neutral' gluon involved. The choice is based then 
on a probability distribution defined by 




where /, J run over all color connections considered. 
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2.4 Interfacing and matching to parton shower programs 

Helac-Phegas generates a Les Houches Accord (LHA) file [33, 34] in a completely au- 
tomatic way, with all the necessary information needed to interface to Pythia and/or 
Herwig parton shower and hadronization programs. 

Possible problems of double counting of jets may arise when interfacing fixed order 
tree level matrix elements to parton shower codes. The reason is that a jet can appear 
both from relatively hard emission during shower evolution or from inclusion of higher 
order matrix elements. In order to deal with this problem a matching or merging algo- 
rithm has to be applied. It removes double counting of jet configurations and provides 
a smooth transition between the part of the phase space covered by parton showers and 
the one described by matrix elements. We refer to [38-43] for a detailed account of these 
algorithms. Here we will briefly present what has been incorporated in Helac-Phegas, 
which is the so-called k± reweighting algorithm based on the k± measure [44,45]. More- 
over, we will also give a brief description of the whole framework, within the so-called 
MLM matching approach [42], which goes beyond Helac-Phegas capabilities and relies 
on the interface with a parton shower algorithm. 

The matching algorithm can be described as follows: 

1. At the beginning one generates all parton level conflgurations for all final state 
parton multiplicities n up to a given N. Due to the presence of collinear and soft 
singularities one has to use a set of parton level kinematical cuts: 

where p^'"'* and are the transverse momentum and pseudorapidity of the final 
state partons, and AR^j is their minimal separation in the {rj, 0) plane. The param- 
eters p™™, rjruax and R^in are called generation parameters, and are the same for all 
n^l, N. 

2. During the generation of events, the renormalization scale (running as) is set ac- 
cording to the CKKW prescription [38,39]. To this end a tree branching structure 
is defined for each event, allowing, however, only for those branchings which are 
consistent with the color structure of the event. In Helac-Phegas this is inherent 
in the matrix element calculation. More precisely for a pair of final-state partons i 
and j, we use the /c^-measure defined by 

k±,ij = ARijmm{pTi,PTj) , 

where ARij = ^ Arjfj + Adfj , while for a pair of initial/final state partons we have 
kj_ = the px of the final state parton. 
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Then if /c^ is minimal the 'partons' are clustered, and the resulting 'parton' is again 
classified as a final state parton with p — Pi +Pj and adjusted fiavor and color fiow. 
In the case when incoming and outgoing partons are clustered, the new 'parton' is 
considered as incoming and its momentum is given by p = pj — Pi- As a result we 
obtain a chain of A;_L-measures. For every node, a factor of Q;s(A;5^ no(jg)/Q;s((5hard) 
multiplied into the weight of the event. For the un-clustered vertices as well as for 
the scale used in the parton density functions, the hard scale of the process Q^ard 
used. No Sudakov reweighting is applied. 

3. At this level, an LHA file has been created, which can be processed independently for 
showering and hadronization. In that sense Helac-Phegas involvement ends here. 
In order to have a self-contained presentation for the reader, we will now give a brief 
description of a possible matching algorithm related to the parton-showered events 
we have already used in Ref. [43]. For all partons resulting from the shower evolution 
a jet cone algorithm is applied. The resulting jets are defined by -EjJmirn ^max 
by a jet cone size i?cius- The parton from the parton-level event is then associated 
to one of the constructed jets. Starting from the parton with the highest pr we 
select the closest jet (1.5 x i?cius) in the pseudo-rapidity/azimuthal-angle space. All 
subsequent partons are matched iteratively to jets. If this is impossible, the event is 
rejected. Additionally, for n < N, matched events with the number of jets greater 
than n are rejected, whereas for n — N, i.e. the highest multiplicity, events with 
extra jets are kept, only if they are softer than the N matched jets. This procedure 
provides the complete inclusive sample. The harder the threshold for the transverse 
energy of the jets used in the matching, E^^, the fewer the events rejected by the 
extra-jet veto (i.e. smaller Sudakov suppression), which means that the role given to 
the shower approximation in the production of jets is more dominant. On the other 
hand using a lower threshold the role of the matrix-element calculation is enhanced, 
which can be seen as equivalent to a more substantial Sudakov suppression, reducing 
thus the role of the parton shower algorithm in the production of the inclusive 
sample. 

In the current version, Helac-Phegas 1.2.0, the fcj^-reweight algorithm described so 
far is included and it can be applied (optionally) to all processes. Moreover, the hard 
scale of the process can be defined by the user in the getqcdscale .h file as described in 
the next section. 
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3 Running the code 



3.1 Description of the code 

The basic structure of the program is already described in [9]. Just to summarize, the 
program reads several input variables and in the so called first phas^ the solution to the 
Dyson-Schwinger equations is constructed. In contrast to other programs all necessary 
information is kept in memory, no extra subdirectories are generated or needecj^. This 
is the reason why Helac-Phegas has a really compact and small source directory. All 
arithmetics in the first phase is completely integer and results in the construction of what 
we call a skeleton, a solution of the recursive equations. In the second phase based on 
the information produced so far, the program generates completely automatically phase 
space points by Phegas, and then calculates the corresponding matrix elements, for each 
color connection configuration separately. In this way not only the total weight can be 
computed but also the weight of each color connection configuration is available for free. 
Color and phase space un-weighting is therefore straightforward. 

Helac-Phegas incorporates all Standard Model particles and couplings in both uni- 
tary and Feynman gauges. Unstable particles are treated in a fully consistent way, either 
by including a fixed width or by using the complex mass scheme [46-48]. Moreover, 
when final states are considered, all correlations (spin, color) are automatically taken into 
account with no approximation involved. 



3.2 Preliminary operations 

Unpacking the code with the command 

tar zxvf helac_l . 2 . . tar . gz 

will create the directory helac_l .2.0. 

The program makes use of the Les Houches Accord PDF Interface library (LHAPDF). 



If this library is not installed on your machine, download it from |http : //projects .hepf orge . org/lhapdf / 
and follow the installation instructionj^. 

After that you have to edit once for all the file myenv. Here you have to set up the 
following environmental variables: 

FC is your FORTRAN compiler; nowadays gf ortran is freely available and can be 

downloaded from http://gcc.gnu.org/wiki/GFortran. The code has been 
also tested with other compilers, like Lahey Fujitsu If 95, Intel Fortran 
if ort and GNU ggfl 



^By default both first and second phases are executed in one run. 

^There is, however, the option to keep this information in a file and use it later on as well. 
^Note that if you do not have root privileges, you can install the library in a local directory. 
^The following compilers are available at IxplusOcern. ch: gfortran, lf95 and g95. 
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FORTRAN_LIBRARIES is the name of the directory in which your FORTRAN libraries are 
located; in some cases they are automatically found by the compiler, and you 
do not need to write them explicitely; 

LHAPDFLIB is the absolute path of your Les Houches Accord PDF Interface library; in a 
standard installation of LHAPDF, it is /usr/local/lib; 

LHAPDFSETS is the absolute path of the directory in which the PDF sets are stored; in a 
standard installation of LHAPDF, it is /usr/local/share/lhapdf /PDFsets. 

In case LHAPDF is not available Helac-Phegas can still run in standalone mode 
(see below the definition of the lhapdf keyword) using the CTEQ6L1 best fit (LO fit with 
LO as). In both cases a running as is available (LHAPDF provides its own running as). 

3.3 How to run a process 

The user interface consists of the following files: 

run.sh is a bash script that reads the input files, compiles Helac-Phegas and runs 
it. run_lsf . sh and run_lsf 2 . sh are its variants to run with the LSF systeir§|. 

default . inp contains a list of keywords, their default values and a short comment. The 
user should not modify this file, but use it just as a reference. 

keywords . list is a list of all the available keywords; it is needed by run.sh. 

user . inp is the only file that the user needs to edit. Here the user can select the process 
and modify the default values of many parameters!^ 

getqcdscale .h this is the file where the definition of the QCD scale to be used in the 
structure function as well as in the as can be set by the user. By default (no 
editing is required) it will use Qhard = Mz- 

sub_lsf and sub_lsf 2 are the bash scripts that send the jobs to the LSF system. 
Let us now illustrate the structure of user, inp with an example. 
# Compulsory information 

colpar 1 # colliding particles: l=pp, 2=ppbar, 3=e+e- 

inist 3-3 # initial state; enter to sum over initial states 

finst 11 -11 # final state 

energy 14000 # collision energy (GeV) 



^http : / /batch . web . cern . ch/batch/ 

^The user is free to rename this file. We will refer to it as user, inp in the following just for clarity. 
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# For reference, here is the particle numbering: 

# ve e u d vm mu c s vt ta t b photon z w+ w- g h f+ f- chi jet 

# 1 2 3 4 5 6 7 8 9 10 11 12 31 32 33 34 35 41 42 43 44 100 

# The respective antiparticles have a minus sign (for example: positron is -2) 

# A jet in the final state is denoted by the number 100 

# Enter here your additional commands if you wish to alter the default values 

With the first four keywords we select the process we want to run and the collision 
energy. In the example, we run uu ti as part of a proton-proton collision at 14 TeV. 

Entering two initial state particles with the inist keyword, we have implicitly chosen 
the single process mode of Helac-Phegas. The other possibility is the summation mode, 
where the initial state particles are hadrons (pp in the example). Helac-Phegas will 
find all the partonic processes that produce the selected final state and will sum over 
them. To achieve this, the user should replace the third line in the example with 

inist # initial state; enter to sum over initial states 

In this example, Helac-Phegas finds 9 partonic subprocesses that produce tt in the 
final state (i.e. the partonic initial state can be qq and qq for q = u,d, s, c, and gg). The 
b quark, by default, is not taken into account as an initial state parton, but this can be 
changed with the qnum keyword which defines the number of quark fiavors considered in 
the initial and final state, as described in the next subsection. 

The summation over partonic subprocesses can also be performed in the final state. 
The particle identification ID=100 which we call a 'jet' has the following meaning: the 
programme will automatically find all parton level subprocesses, replacing this ID=100 
with either a gluon or a quark (anti-quark). Which quark flavors are taken into account 
in this subprocesses generation is again defined by the qnum keyword. For example, if we 
replace third and fourth line in the example with 

inist # initial state; enter to sum over initial states 

finst 11 -11 100 # final state 

Helac-Phegas will look for all partonic subprocesses having ti + 1 parton in the final 
state. By default, the t quark is not taken into account as a 'jet' initiator. 

In summation mode, the programme calls two auxiliary scripts, script, sh and 
script, l.sh, located in the Summation_Processes subdirectory. 

In case of many jets in the final state, the number of subprocesses increases rapidly. 
Since the generation for different subprocesses is independent, the program is trivially 
parallelizable. As an example, we have implemented suitable scripts using the LSF system 
at lxplus@cern.ch, see subsection 13.51 for more details. 
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There are many parameters and options that can be changed by the user. All of them 
are listed in the default, inp file and described in more detail in subsection 13. 4[ 

For example, the line of default . inp referring to the number of Monte Carlo iterations 

is 

nmc 100000 # number of montecarlo iterations (single mode) 

To change the number of Monte Carlo iterations from 100000 to let say 500000 one has 
to add the line 

nmc 500000 

at the end of the user . inp file. Choices made here will overwrite the default values. 

After editing the user . inp file, you can compile and run the program by typing 
./run.sh user. inp 
or 

./run.sh user. inp myenv-xxx. 

In this last case, a file myenv-xxx will be use instead of myenv. 

Summarizing, the user can edit user, inp but is not allowed to change its structure. 
In the upper part of the file (lines 1-13) the only thing the user can do is changing the 
numerical values of the first four keywords to select the process. In the lower part of the 
file (starting from line 14) the user can 

1. Add lines in the form keyword value to alter the default values; 

2. Add empty lines or comment lines. 

In case of wrong keywords, repeated keywords or invalid values, the script will complain 
and ask the user if she/he wants to continue anyway. 

3.4 Keywords explanation 
3.4.1 Compulsory keywords 

The following keywords must be present in each user, inp file (lines 2-5) in order to 
select a process which should be run. The symbols, s, n, b and x correspond to the string, 
integer, boolean (logical) and real data types respectively. 

• colpar n: Colliding particles. 

Set to n = 1 for a pp collision, n = 2 for pp, n = 3 for e"^e^. 

• inist nin2: Particles in the initial state. 

This keyword implicitly selects the running mode. If both rii and n2 are chosen 
from Table [H Helac-Phegas will run in single process mode. If rii = and n2 is 
omitted, Helac-Phegas will run in summation mode. 
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• f inst ni n2 . . . nfc : Particles in the final state. 

Allowed values are the numbers listed in Table [H In summation mode, = 100 
represents a jet. 

• energy x : Collision energy in GeV. 
3.4.2 Optional keywords 

The following keywords allow the user to change many default parameters and options. 
They can be entered beginning from the 14th line of user.inp. 

General options 

• exec s: Name of the executable file to be produced. 
Default name is helac-phegas.exe. 

• outdir s : Name of the output directory. 
Default name is RESULTS. 

• nglu n: Number of gluons in the final state. 

In summation mode and with jets present in the final state, select only the subpro- 
cesses having n gluons in the final state. Set to n = — 1 to not constrain the number 
of gluons in the final state. 
Default value is —1. 

• qnum ni?T,2: Number of quark flavors in the initial and final state. 

In summation mode, when looking for the allowed subprocesses Helac-Phegas 
will take into account only the ni lightest flavors in the initial state, and the n2 
lightest flavors as jet initiators in the final state. 
Default values are ni = 4 (u, d, c, s) and n2 = 5 {u, d, c, s, b). 

• ktrw b: Switch for the k± reweighting algorithm. 

See section [2^41 for more details. Set F, f or to switch off the algorithm, or set T, 
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Number 
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11 
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32 


Z 


44 


X 



Table 1: Particle numbering 
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t or 1 to switch it on. 
Default value is false. 



• Irgnc b: Switch for large Nc limit. 

In the limit Nc ^ oo, only the diagonal terms of the color matrix Cjj, cf. Eq. ( fTTll . 
survive, and all O^N'"^) terms are neglected. The interferences between different 
colour flows vanish which speeds up the calculation. This option can be used e.g. 
to test, whether the large limit approximation can reliably describe the physical 
procesJil. Moreover, in the case when ktrw is switched on, the user still has the 
option to switch on or off the large Nc limit. 
Default value is false. 

• error s: Name of the file, where, every 1000 iterations, the cross section and its 
Monte Carlo error are recorded. 

Default is err_f ile. 

• gener n: Choice of the phase space generator. 



Besides the default Phegas [15] described in Section ETT], there are two other possi- 
bilities. The first one is Rambo [49], a fiat phase-space generator. It generates the 
momenta distributed uniformly in phase space so that a large number of events is 
needed to integrate the integrands to acceptable precision, which results in a rather 
low computational efficiency. The second one is Durham, an algorithm that for a 
given order of the final state particles expresses the phase space in terms of the /c^, 
and y variables, where y is the difference of the rapidities of two nearest particles. 
Adaptation over all possible orderings of the final state particles is also included. 
For more details see Ref. [12]. Set to n = to use Phegas, n = 1 for Rambo or 
n = 2 for Durham. 
Default value is 0. 

• repeat n: Helac-Phegas run option. 

Set to n = to run both phases of Helac-Phegas, n = 1 to run only the first 
phase, where the skeleton for a given process is constructed, or n = 2 to run only 
the second phase. 
Default value is 0. 

• ranhel n: Choice of the set-up for the helicity configurations. 

There are two possibilities, either exact summation over all helicity configurations 
is performed or Monte Carlo integration is applied [9]. For example for a massive 
gauge boson the second option is achieved by introducing the polarization vector 



''^Let us stress here, that default calculation in Helac-Phegas takes into account the full color struc- 
ture, no approximation is involved. 



(14) 
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where e (0, 2%) is a random number. By integrating over (f) the correct sum over 
helicities is obtained 

A — i 

The same idea is applied to the helicity of (anti)fermions. Set to n = to calculate 
explicitly all helicity configurations or n = 1 to do a Monte Carlo over helicities. 
Default value is 1. 

alphasrun n: Switch for the running of cc^. 

Set to n = to have a fixed value oi ag oi n — 1 to have a running coupling 
OisiQ^)- In the later case the Q scale has to be implemented by the user in the 
getqcdscale .h file, where a few examples have already been provided. 
Default value is 0. 

gauge n: Choice of the gauge. 

Set to n = for the Feynman gauge or n = 1 for the unitary gauge. 
Default value is 1. 

ihiggs b: Switch for the inclusion of the contribution coming from the Higgs boson. 
Default value is false. 

widsch n: Scheme for the introduction of the width of W, Z bosons as well as the 
top quark [46-48]. Set to n = for the fixed width scheme or n = 1 for the complex 
mass scheme. 
Default value is 0. 

qcd n: Option for the type of interactions. 

Set to n = to have only electroweak interactions, n = 1 to have both electroweak 
and QCD and n = 2 for QCD only. 
Default value is 1. 

unwgt b: Switch for the un-weighting procedure. 

If set to false events are not unweighted. If set to true events are unweighted. 
Default value is true. 

preunw n: Pre-unweighted events. 

Number of events generated to calculate the maximal weight to start the un- 
weighting procedure. It is recommended that the user tries different values in order 
to obtain a better un-weighting efficiency. 
Default value is 10000. 

unwevt n: Number of unweighted events to be generated. 

In the summation mode, n is the total number of unweighted events, which is 
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distributed between all possible parton level processes according to their relative 
contribution to the total cross section. 
Default value is 50000. 

nmc n: Number of Monte Carlo iterations in the single process mode. 
Not used in summation mode. 
Default value is 100000. 

nmcl n: Number of Monte Carlo iterations in the first step of the summation mode. 
Not used in single process mode. 
Default value is 100000. 

nmc2 n: Number of Monte Carlo iterations in the second step of the summation 
mode. Not used in single process mode or if un-weighting is switched off. 
Default value is 3000000. 

optim 711^2 713^4 715^6: Optimization options for the Phegas phase space gen- 
erator. Optimization is performed when the program reaches the Monte Carlo 
iterations 

{oi = Til, 02 = 0i + n2n3, 03 = 02 + n2nl, ...,0^^ = o^^-i + narig^"^} . 

It stops when the maximum number of optimizations is reached, or when the 
7i5-th Monte Carlo iteration is reached. Uq is a flag: if set to optimisation is not 
performed, if set to 1 it is performed. 
Default values are 10000 10000 1 8 100000 1. 

lhapdf h: Switch for the use of the Les Houches Accord PDF interface library. 
If set to F, Helac-Phegas will run in standalone mode, using the CTEQ6L1 best 
fit (LO fit with LO Us). 
Default value is false. 

pdf n: Choice of the PDF set. 

Used only if lhapdf is set to true, n must be equal to the LHAGLUE number of 
the desired PDF set (see PDFsets. index or Ref. [28]). Enter 7i = to run without 
PDFs. 

Default value is 10042, which corresponds to CTEQ6L1 (LO fit with LO as). 

pythia 7ii7i2: IDWTUP and NPRUP parameters for the LHA sample files. 
IDWTUP dictates how the event weights are interpreted by the showering and hadroniza- 
tion event generators. NPRUP is the number of different user processes, see [33] for 
more details. 

Default values are 7ii = 3, which corresponds to an option with unweighted events 
with the weight = +1, and 7i2 = 1, which corresponds to the one process only. 
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• constants b: if 6 = 1 the user can provide its own constants. h file in accordance 
to the prototype given in the constants_std.h file. 
Default value is 6 = 0, see below the description of physical constants. 

pp(bar) cut^ 

The following options are read only if colpar is set to 1 or 2 (respectively, pp and pp). 
All dimensional quantities are in GeV. When the word 'quark' appears, light quarks and 
gluons have to be understood. 



Key woici 




ClcidUll 


cutof fp X 


J- XT 1 

Luton value 


l.Ud-3 


minptl X 


Minimum lepton pt 


l.Ud-o 


minptq x 


Minimum quark Pt 


8.0 


minptb x 


Minimum bottom Pt 


8.0 


minptt X 


1\ /T * * J- 

Minimum top Pt 


0.0 


minptp X 


1\ /T * * 1, J- 

Minimum photon pt 


8.0 


maxrapl x 


Maximum lepton rapidity 


10.0 


maxrapq x 


Maximum quark rapidity 


2.0 


maxrapb x 


Maximum bottom rapidity 


2.0 


maxrapt x 


Maximum top rapidity 


10.0 


maxrapp x 


Maximum photon rapidity 


2.0 


mindrll x 


Minimum AR separation between lepton and lepton 


0.0 


mindrlq x 


Minimum AR separation between lepton and quark 


0.0 


mindrqq x 


Minimum AR separation between quark and quark 


0.7 


mindrqb x 


Minimum AR separation between quark and bottom 


0.7 


mindrbb x 


Minimum AR separation between bottom and bottom 


0.7 


mindrpf x 


Minimum AR separation between photon and fermion 


0.0 


minmqqp x 


Minimum quark-quark invariant mass in the pp/pp case 


0.0 


minmqb x 


Minimum quark-bottom invariant mass 


0.0 


minmbb x 


Minimum bottom-bottom invariant mass 


0.0 



e+e" cutqil 

The following options are read only if colpar is set to 3 (e+e~). All dimensional quantities 
are in GeV. 

^Cuts not already provided may be added in the 2.cuts_auto_pp.f file. 
^Cuts not already provided may be added in the 2.cuts_auto_ee.f file. 
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keyword 


description 


default 


cutof f e X 


Cutoff value 


l.Od-3 


minenl x 


Minimum lepton energy 


10.0 


minenq x 


Minimum quark energy 


10.0 


minenp x 


Minimum plioton energy 


10.0 


minanglb x 


Minimum angle (degrees) between lepton and beam 


5.0 


minangqb x 


Minimum angle (degrees) between quark and beam 


5.0 


minangpb x 


Minimum angle (degrees) between photon and beam 


5.0 


minangll x 


Minimum angle (degrees) between lepton and lepton 


5.0 


minanglq x 


Minimum angle (degrees) between lepton and quark 


5.0 


minangqq x 


Minimum angle (degrees) between quark and quark 


5.0 


minangpf x 


Minimum angle (degrees) between photon and fermion 


5.0 


minmqq x 


Minimum quark-quark invariant mass in the e~^e~ case 


10.0 



Physical constants 

The user can set the values of many physical constants using the keywords reported in 
Table El 

By default the electroweak couplings are given by 

y 9weak = {4:V2GFy^'^mw, Oiem = V2GF^T^^Sm^6'^„/7^ 

niz ) 

However, the user can optionally give values for aem (alphaem) and siri^dw (sin2thetaw). 

For more advanced users, it is recommended to edit and provide their own constants . h 
file. A prototype is included, see the constants_std.h file. The CKM [50,51] matrix is 
also possible to be defined in this file, as well as masses and widths for all particles. 

3.5 How to run HELAC using the LSF system 

Running Helac-Phegas on a single machine can become impractical when the final 
state contains many particles and jets, either because of time constraints or because of 
the amount of output that can be produced. A possible solution to this problem is to run 
the n involved subprocesses at the same time on a computer cluster. We provide scripts 
to run Helac-Phegas on the CERN's batch system 

The keywords queue 1 and queue2 allow the user to choose the queue for the subpro- 
cesses in the first and second phase of the Helac-Phegas run. 

By default, jobs are sent to the 8nm queue, 8 minutes, which is suitable only for short 
runs. Two options are available: 

1. Write everything, also the intermediate files, which can be quite large, in the AFS 
area. 

In this case, a user . inp file has to be prepared as usual, and then 
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keyword 


description 


default 


gfermi x 


Fermi coupling constant 


1.16639d-5 


sin2thetaw x 


Sinus squared of Weinberg angle siri^Ow 


-1* 


alphaem x 


Electromagnetic coupling constant 


-1* 


alphas2 x 


Fixed strong coupling constant, used only if alphasrun is 


0.118 


zmass X 


Z mass 


91.188 


zwidth X 


Z width 


2.446 


wmass X 


W mass 


80.419 


wwidth X 


W width 


2.048 


higmass x 


Higgs boson mass 


130.0 


higwidth X 


Higgs boson width 


4.291d-3 


emass x 


e mass 


0.0 


mumass x 


mass 


0.0 


taumass x 


T mass 


0.0 


umass X 


u quark mass 


0.0 


dmass x 


d quark mass 


0.0 


smass X 


s quark mass 


0.0 


cmass X 


c quark mass 


0.0 


bmass x 


h quark mass 


0.0 


tmass X 


t quark mass 


174.3 


twidth X 


t quark width 


1.6 



Table 2: Keywords and default values regarding some physical constants. 

*The value —1 for the keywords sin2thetaw and alphaem means that they are defined 

as described in the text. 

. /sub_lsf 

should be run. The user will be asked to provide some information like job name, 
queue, user input file and environment file. The queue will only be used for the 
controller job, and may, therefore, be short. 

2. Write everything on the remote, Ixb, hosts (suggested method). 

At the end of the run, the sample will be written in the AFS directory from which the 
job has been submitted. The usual output directory will be written, in compressed 
form, in the user's personal CASTOR area. 

Similarly to the previous case, a user . inp file has to be prepared and then 
./sub_lsf2 

should be run. The user has to provide a job name, queue, user input file and 
environment file. The script will automatically resubmit a job to a longer queue, if 
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it failed. Notice, however, that this is only useful in the case of a failure due to time 
constraints. 

There are two reasons for using this method: 

1. To save time. More subprocesses run at the same time, while with the standard 
procedure they run one after the other. 

2. To save space. The output files produced can exceed your disk quota. 

The first procedure (sub_lsf) saves time but not space, whereas the second one 
(sub_lsf 2) saves also space. Another reason to prefer the second procedure is the auto- 
matic resubmission of failed jobs, which only is implemented in this case. 

3.6 Description of the output 

All the outputs produced by the Helac-Phegas program will be written in a subdi- 
rectory, whose name is defined by the keyword outdir. Let us distinguish two main 
cases. 

3.6.1 Single process mode 

In the output directory you will find the following files. 

• Executable files. The Helac-Phegas executable (whose name is defined by the 
keyword exec) and, only if unwgt is true, the executable for the un-weighting pro- 
cedure, unwei.exe. 

• User input file, called user.inp throughout this paper. 

• input_sp and output. Respectively input and output redirected files of the Helac- 
Phegas executable. If unwgt is true, the output of unwei.exe is appended to 
output. 

• The error file, whose name is defined by the keyword error. Every 1000 iterations, 
the cross section and its error are recorded to this file. The three columns are 
respectively the number of iterations, the cross section and the error. 

• The kinematic file, kine_XXX . out, where XXX represents the process. Used for 
debugging. 

• The event sample file, sampleXXX . Ihe, where XXX represents the process, in the 
standard format dictated by the Les Houches accord [33,34]. It is generated only if 
unwgt is set to true. 
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3.6.2 Summation mode 



In the output directory you will find 

• The Helac-Phegas executable. 

• User input file. 

• infile, the input file for processes. f, that generates all the subprocesses. 

• The (gzipped) sample file, sampleXXX . Ihe . gz, where XXX describes the final state. 
It is generated only if unwgt is set to true. 

• A subdirectory 1, where all the outputs generated in the first phase, generation of 
weighted events, will be written. 

• A subdirectory 2, where all the outputs generated in the second phase, generation 
of unweighted events, will be written. It is generated only if unwgt is set to true. 

The subdirectory 1 contains 

• input, a skeleton of input files, needed to generate the real input files for each single 
subprocess. 

• output, the output of readoutput . f . 

• results . output. Each group of two lines refers to a single subprocess. The first 
line is the cross section and the second is the error. 

• cross_sections. The n-th line shows the percentage of the total cross section 
represented by the n-th subprocess. 

• Subdirectories directory _n, one for each subprocess. The content of each of these 
directories resembles the one described for the single process mode, except that most 
of the files are gzipped to save space. 

The content of the subdirectory 2 is analogous to 1, with the main exception that in 
the results . output file a third and a fourth line for each subprocesses is added, showing 
information on generated unweighted events. 



3.6.3 Output of the main Helac-Phegas executable 

A substantial part of the information contained in the output of helac-phegas . exe is self- 
explainable. In this section we would like to make the user familiar with the main points 
of this file. Lets take as the basic example the following process u(pi)u(p2) — > gipajgip^}- 
The output file begins with the reproduction of the input file, followed by masses and 
widths of all SM particles. Then we have 
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QCD INCLUDED, g= 1.217715784776720 0.1180000000000000 
UqUaGOGO 

TIME= 2007 10 8 180 19 26 52 966 
avhel , avcol , symet 

0.2500000000000000 0.1111111111111111 0.5000000000000000 

where the process is clearly stated UqUaGOGO, together with the QCD coupling and factors 
that correspond to the average over initial helicities, colors and the symmetry factor. 
Then the 'first phase' starts. For each color connection configuration the color connection 
assignments for all particles is given, namely 

the colour of particles ONE 2 3 1 
the colour of particles TWO 12 3 

the first line is the color connection the second the anti-color one (kept fixed at the ordinal 
permutation). Then the solution for the Dyson-Schwinger equations is presented in the 
following form 
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It simply encodes all information - to be used internally - corresponding to the fusion 
of an u with momentum p2 with the g with momentum ps to produce an u with total 
momentum p2 +P3 (all momenta taken incoming). At the end of this part a line with 

the number of Feynman graphs = 3 
the number of Feynman graphs = 3 

followed by the color matrix is given. The number of Feynman graphs for this process is 
of course three. Then the number of channels to be used by Phegas (always the number 
of graphs +1) is given with the number of MC points nmc and optimization parameters. 

Number of channels 4 
Number of MC points 



1 
1 
3 
3 
1 
1 
3 
3 
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nmc = 100000 
nopt ,nopt_step, optf ,maxopt , iopt 
100 10000 1.000000000000000 8 1 
number of channels= 4 
NUMBER OF CHANNELS 4 

Since Phegas is using Feynman graphs as the basic structure of the multichannel phase- 
space mappings, a description of each Feynman graph is given as follows: 

the graph 1 

14 -3 12 35 2 -3 

12 35 4 35 8 35 

After a few lines with self-explained information the cuts 



the cuts 



pt of 3 




particle 


8, 


.000000000000000 


energy of 3 




particle 


8, 


.000000000000000 


pt of 4 




particle 


8, 


.000000000000000 


energy of 4 




particle 


8, 


,000000000000000 


rapidity of 


3 


particle 




2.000000000000000 


rapidity of 


4 


particle 




2.000000000000000 


cos-beaml of 


3 


particle 




0.9640275800758169 


cos-beaml of 


4 


particle 




0.9640275800758169 


cos-beam2 of 


3 


particle 




0.9640275800758169 


cos-beam2 of 


4 


particle 




0.9640275800758169 


DR 3 with 


4 0.7000000000000000 


cos of 3 with 


4 0.7648421872844885 


mass of 3 with 4 5.486364919287221 



applied are given for each particle or pair of particles. Then a bunch of lines follows in 
the form 

sigma= 0.175490D+03 0.709135D-01 491 1000 1000 



where the calculated cross section sigma with the percentage error, the number of events 
passing the cuts, the number of phase-space points used as well as the number of phase- 
space points tried, is given. 

At the end of the file one normally has 
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out of 100000 100001 points have been used 
and 50419 points resulted to =/= weight 
whereas 49582 points to weight 

estimator x: 0.182757D+03 

estimator y: 0.147435D+01 

estimator z: 0.124416D-02 

average estimate : 0.182757D+03 
+\- 0.121423D+01 

variance estimate: 0.147435D+01 
+\- 0.352726D-01 
Iwri: points have used . OOOOOOOOOOOOOOOE+00 
2212 2212 7000.000000000000 7000.000000000000 3 1 
I error: 0.6643944554436630 

which is translated to the total cross section in nb . 182757D+03 the MC error also 
in nb 0.121423D+01 and the percent error 0.6643944554436630%. 

If the unwgt is set to true the distribution of the weights is plotted for the number 
of events given by preunw. Then the un-weighting procedure starts and the number of 
unweighted events is additionally printed. In case of un-weighting the final histogram of 
the weights is also printed. 



3.7 External phase space points 

In this mode the user can provide by herself/himself momenta for the particles participat- 
ing in the process under consideration and get back information on the squared matrix 
element summed over spin and color degrees of freedom. Additionally, the color ordered 
amplitude for a given color connection and helicity configuration can be obtained. 

• onep : Set to true in order to have the described option. 
Default value is false. 

• mom : Name of the input file provided by the user with the momenta of the particles, 
in the format Px,Py,Pz, E,m,w, where w is a weight for phase space generation, 
usually set w — 1. A prototype is included under the name mom. 

Default name is mom. 

• momout: Name of the output file. 
Default name is momout. 



3.8 Benchmarks 

From the physics point of view Helac-Phegas has already been used in several contexts 
to produce results in physics. A few examples are: the Monte Carlo generator Nextcal- 
IBUR [52] that includes the ExCALlBUR [53] phase space generation and Helac-Phegas 
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matrix elements, used in e^e~ analysis; the extensive comparison with Sherpa on six- 
fermion production processes [54]; and more recently the participation in the W + n 
jets comparative study [43] of different Monte Carlo codes namely Alpgen, Ariadne, 
Helac-Phegas, MadEvent and Sherpa. 

On the technical side, we provide a number of benchmark calculations in the web-page 
of the code, so that the potential user can test and validate his own results. 

4 Outlook 

The current version is named 1.2.0; it is meant to be publicly used: in that sense we wel- 
come any bug report or simple question. The first number introduces major changes, the 
first one described below. The second number refers to minor changes and/or additions, 
whereas the last one to bug corrections. 

A version with all Higgs-gluon and Higgs-photon couplings, in the large rritop limit will 
follow as 1.3.0. It is already tested and will be available soon on the Helac-Phegas 
web-page. 

A new version suitable for processes with many colored particles, i.e. with a number 
of equivalent gluons more than 9, for instance gg —>■ 8g and more, has been already 
developed and tested [12]. It will be incorporated in version 2. 

We are also working on the inclusion of MSSM particles and couplings. 
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